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We develop equilibrium and kinetic theories that describe the assembly of viral capsid proteins on a charged 
central core, as seen in recent experiments in which brome mosaic virus (BMV) capsids assemble around 
nanoparticles functionalized with polyelectrolyte. We model interactions between capsid proteins and nanopar- 
ticle surfaces as the interaction of polyelectrolyte brushes with opposite charge, using the nonlinear Poisson 
Boltzmann equation. The models predict that there is a threshold density of functionalized charge, above which 
capsids efficiently assemble around nanoparticles, and that light scatter intensity increases rapidly at early times, 
without the lag phase characteristic of empty capsid assembly. These predictions are consistent with, and en- 
able interpretation of, preliminary experimental data. However, the models predict a stronger dependence of 
nanoparticle incorporation efficiency on functionalized charge density than measured in experiments, and do 
not completely capture a logarithmic growth phase seen in experimental light scatter. These discrepancies may 
suggest the presence of metastable disordered states in the experimental system. In addition to discussing fu- 
ture experiments for nanoparticle-capsid systems, we discuss broader implications for understanding assembly 
around charged cores such as nucleic acids. 



I. INTRODUCTION 

The cooperative assembly of heterogenius building blocks 
into ordered structures is crucial for many biological pro- 
cesses, such as the assembly of protein subunits around a vi- 
ral nucleic acid to form a capsid. Understanding how nucleic 
acid-protein interactions drive cooperative assembly could en- 
able antiviral drugs that block this essential step in viral repli- 
cation. At the same time, engineered structures in which 
viral capsids assemble around synthetic cargoes show great 
promise as delivery vehicles for drugs or imaging 

agents flIslSfl and as subunits or templates for the synthesis 
of advanced multicomponent nanomaterials flSllRllll 
Realizing these goals, however, requires understanding at a 
fundamental level how properties of cargoes and capsid pro- 
teins control assembly rates and mechanisms. With that goal 
in mind, this article presents coarse-grained thermodynamic 
and kinetic models that describe experiments in which brome 
mosaic virus (BMV) capsid proteins assemble around spheri- 
cal charge-functionalized nanoparticle cores. 

Identifying mechanisms for simultaneous assembly and 
cargo-encapsidation. Elegant experiments have studied the 
assembly of capsids around central cores consisting of nu- 
cleic acids (e.g. Jisjl^jisl [jil O [ijl [ij, |2o|], inorganic 




2311 . and charge-functionalized cores 
27|,|28|]. However, assembly mechanisms of 
viruses and virus-like particles remain incompletely under- 
stood because cargo-protein interactions and the structures of 
transient assembly intermediates are not accessible to experi- 
ments. For the case of empty-capsid assembly, Zlotnick and 



coworkers achieved great insight about assembly mechanisms 

with simplified theories (e.g. H S H H H H H), 

for which a small number of parameters are fit by compari- 
son with experimental assembly data over a range of protein 
concentrations and pH values. The theories developed in this 
work enable similar comparisons for the assembly of capsids 
around electrostatic cores, and could thereby elucidate mech- 
anisms for simultaneous assembly and cargo encapsidation. 

Prior theoretical works have led to important insights about 
assembly around polymeric cores, but were equilibrium stud- 



ies 1136, 
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40l ? ] or considered specific RNA- 



mediated assembly pathways with phenomenological descrip- 
tions of protein-nucleic acid interactions 14 ll 14211 . Computer 
simulations of a model in which subunits assemble around 



rigid spherical cores ll43ll44ll suggest that a core with a geom- 
etry commensurate with the lowest free energy empty-capsid 
morphology can increase assembly rates and the e fficiency of 
assembly (compared to empty capsid assemblyll45 
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541]) through the effects of heterogeneous 
nucleation and templating; however, assembly is frustrated for 
core-subunit interaction strengths that are too strong|j43, 551. 
Cores that are not size-matched with the lowest free energy 
capsid morphology can direct assembly into alternative capsid 
morphologies for optimal core-subunit interaction strengths 
1 4311 . Important limitations of the simulations are that ex- 
tensive comparison with experimental data is computationally 
intractable, and the simulated models do not explicitly rep- 
resent electrostatics or the polymeric nature of disordered N- 
terminal tails in capsid proteins. 

The present work aims to overcome these limitations 



2 



by calculating the electrostatic interaction between charge- 
functionalized nanoparticles and capsid proteins, and then 
constructing simplified thermodynamic and kinetic theories 
that describe the simultaneous assembly and cargo encapsi- 
dation process. We explore the effect of cargo-capsid protein 
interactions on assembly mechanisms by presenting predicted 
assembly rates and nanoparticle incorporation efficiencies as 
functions of the functionalized charge density, nanoparticle- 
protein stoichiometric ratio, and the capsid protein concentra- 
tion. Because each of these parameters can be experimentally 
controlled, the theoretical predictions suggest a series of ex- 
periments that can elucidate the relationship between cargo- 
capsid protein interactions and assembly mechanisms. 

Although there is currently not sufficient experimental data 
to estimate values for several unknown parameters (adsorption 
and assembly rate constants, and subunit-subunit binding en- 
ergies), we qualitatively compare the theoretical predictions to 
preliminary experimental time-resolved light scattering data 
(Fig. [T^) and measurements of the efficiency of nanoparticle 
incorporation in well-formed capsids (Fig. (TJi). With physi- 
cally reasonable values for the unknown parameters, the pre- 
dicted light scatter increases rapidly at short times without a 
lag phase, consistent with the experimental measurements, but 
the theory does not completely capture the functional form of 
light scatter at long times. The theory predicts that nanoparti- 
cles are incorporated into capsids only above a threshold func- 
tionalized charge density, which is similar to the lowest charge 
density for which significant incorporation is observed in ex- 
periments. However, predicted incorporation efficiencies rise 
to 100% over a narrow range of charge densities, whereas the 
experimental data shows a gradual increase in the incorpo- 
ration efficiency with increasing charge density. We discuss 
experiments and potential improvements to the theory which 
might elucidate these discrepancies. 

This paper is organized as follows. In Section|lI]we calcu- 
late the interaction between capsid proteins and functionalized 
charge on nanoparticle surfaces. In Section Hill we present 
a model for the thermodynamics of capsid assembly in the 
presence of rigid spherical cores, followed in Section HVlbv 
a kinetic theory that describes the simultaneous assembly of 
capsid subunits on nanoparticle surfaces and in bulk solution. 
We analyze kinetics and the assembly pathways predicted by 
the model in section |Vl and discuss the connection to current 
and future experiments in Section[yT] We also compare theory 
predictions to simulation results in Appendix A. 
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FIG. 1: (a) Time-resolved light scatter for capsid assembly on 
nanoparticles with a functionalized surface density of ctc = 3 acid 
groups/nm^. (b) Packaging efficiency, or the fraction of nanopar- 
ticles incorporated into T=3 capsids (measured from TEM micro- 
gaphs), varies with ac- The capsid protein-nanoparticle stoichiomet- 
ric ratio is rcp:NP = 1.33, core diameters are 12 nm (which pro- 
mote T=3 capsid formation), the capsid protein concentration is 0.4 
mg/ml (Cp = 10 fiM dimer subunit), and there is 100 mM 1:1 salt 
and 5 mM 2: 1 salt. Packaging efficiencies were measured ait — W 
seconds. The surface density of acid groups is estimated from the 
fraction of carboxylated TEG molecules, assuming a total surface 
density of 3 TEG molecules/nm^ on the nanoparticle surface. Data 
thanks to Bogdan Dragnea (5^ . 



Experimental motivation. This work is motivated by ex- 
periments in which BMV capsid proteins assemble around 
gold nanoparticle cores 1I24' 



2711 functionalized with thio- 



lalkylated tetraethylene glycol (TEG) molecules, a fraction of 
which are terminated with carboxylate groups ll27ll . With 12 
nm nanoparticles the TEG-nanoparticle complexes have a di- 
ameter of roughly 16 nm, which is nearly commensurate with 
the interior of a T=3 BMV capsid (17-18 nm [57]). Ionized 
carboxyl groups provide a driving force for capsid proteins to 
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FIG. 2: The model system. Surface functionalization molecules 
(TEG) are end-grafted to an impenetrable gold sphere. The cationic 
capsid protein N-terminal arms are modeled as polyelectrolytes 
grafted to the inner surface of a sphere (the capsid), which is im- 
permeable to TEG and polyelectrolyte but permeable to solvent and 
ions. 



adsorb and assemble on nanoparticle surfaces. Time-resolved 
light scatter measurements (Fig. [T^), which monitor assem- 
bly kinetics, are strikingly different than light scatter of empty 
capsid assembly at short times (e.g. ifsS, 59, 60]); instead of 
a lag phase (~ 10 seconds - minutes) there is a rapid (^ 1 
second) rise in signal. As described in Section|Vl the theoret- 
ical predictions in this work begin to explain this observation. 
Assembly effectiveness is monitored by using TEM to mea- 
sure packaging efficiencies, or the fraction of nanoparticles in- 
corporated within well-formed capsids. As shown in Fig.[TJ), 
packaging efficiencies increase from 0% to nearly 85% as the 
fraction of TEG molecules that are carboxylate-terminated in- 
creases from 20% to 100%. 

While the light scattering data was obtained from an as- 
sembly reaction entirely at pH = 5, the packaging efficiency 
data shown in Fig. [T^ were obtained with a protocol in which 
capsid subunits were first dialyzed against a buffer at pH = 
7.4, followed by buffer at pH = 5. In this work we con- 
sider only pH = 5, since at pH = 7.4 disordered protein- 
nanoparticle complexes form but well-formed capsids are not 
observed without subsequently lowering pH. Assembly may 
not be favorable at pH ~ 7.4 because binding interactions be- 
tween capsid proteins are less favorable than at lower pH IssI : 
also strong TEG-protein interactions might render disordered 
states kinetically or thermodynamically favorable. 



II. SUBUNIT ADSORPTION ON NANOPARTICLE 
SURFACES 



For many capsid proteins, there is a high density of basic 
residues on flexible N-terminal tails, with numbers of charged 
amino acids that range from six to 30 ll36ll . These positive 
charges help drive capsid assembly in vivo by interacting with 




FIG. 3: The volume fractions of (a) polyelectrolyte segments and (b) 
carboxylate groups in the grafted TEG layer are shown as a func- 
tion of layer (Az = z ~ Rc/l + 1) above the nanoparticle surface 
for several surface acid group densities, with the nanoparticle radius 
Rc ~ 6 nm and the layer dimension I — 0.35 nm. 



negative charges on nucleic acids, and the interaction ofpoly- 
meric tails and nucleic acids is explored in Refs. 136, 42, 611. 
In this section, we analyze the interactions between capsid 
protein tails and surface functionaUzed TEG molecules. For 
simplicity, we neglect other effects, such as hydrophobic in- 
teractions, that might drive subunit adsorption, and we assume 
that all charges on capsid proteins are located in the polymeric 
tails. For the remainder of this work, we will reserve the word 
'polyelectrolyte' to apply to the charged polymeric tails on 
capsid proteins; we will not use this word to refer to the func- 
tionalized TEG molecules. 

Model system. We first consider a dilute solution of cap- 
sid protein subunits with density pp, and a single nanoparticle; 
we consider a finite concentration of cores in Section |lll] At 
high surface coverage or in an assembled capsid the folded 
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portion of capsid proteins acts as an impenetrable barrier to 
the polymeric tails. We therefore represent the nanoparticle- 
capsid complex as a concave polyelectrolyte brush (the inte- 
rior surface of the capsid) interacting with a convex brush of 
the opposite charge (the TEG layer); the model system geom- 
etry is illustrated in Fig. |2] We consider surface functional- 
ization molecules grafted to the outer surface of an impene- 
trable sphere with radius 6 nm. The TEG surface function- 
alization molecules used in Ref. |24| consist of a 10 carbon 
alkanethiol chemically linked to a polyethylene glycol chain 
of 4 monomer units. The surface density of functionalized 
charge groups, ac, is experimentally controlled by terminat- 
ing a fraction TEG molecules with carboxylates, with the re- 
mainder terminated with neutral groups po). We represent 
TEG molecules as freely jointed chains with 9 neutral seg- 
ments end-grafted to the surface and one charged segment 
(with variable valence) at the free end. 

The N-terminal tails, or 'arms' of BMV capsid proteins 
contain nine positive charges in the first 20 amino acids and 
the first 25-40 residues are disordered in capsid crystal struc- 
tures [15711 . Since BMV capsids assemble from protein dimer 
subunits 1621 l63l 16411 . in our model each subunit has a net 



charge of = 18. We represent the two N-terminal arms 
from each dimer subunit as a single polymer of 36 segments 
with Kuhn length 3.5 A, and valence 0.5. These polyelec- 
trolytes are end-grafted to the inner surface of a sphere with 
radius 9 nm, the inner diameter of a T=3 BMV capsid ll57ll . 
The model capsid is permeable to solvent and ions, but im- 
penetrable to TEG and polyelectrolytes segments. The lat- 
ter condition is reasonable for high surface coverages of cap- 
sid subunits and in the strong coupling limit, in which case 
the polyelectrolyte configurational entropy is negligible com- 
pared to electrostatic effects. For low charge densities, the 
results do not change qualitatively if the impermeable capsid 
is eliminated, but the equilibrium concentration of adsorbed 
polyelectrolytes increases by up to 50%. 

Numerical calculations. To our knowledge, interacting 
polyelectrolyte brushes with this geometry have not been stud- 
ied, although the interaction of capsid arms with a nucleic 
acid is considered in Refs. 136116511 . For many cases we con- 
sider, the electrostatic potential in the region of the capsid and 
nanoparticle surface is large; therefore, we will solve the non- 
linear Poisson Boltzmann (PB) equation. To do so, we employ 
the method of Scheutjens and Fleer (SF) I66ll67ll68ll . in which 



ter, as well as the dissociation equilibrium for carboxylate and 
water groups, are solved numerically with a self consistent 
field approximation on a lattice. The calculation accounts for 
the finite size of all species and the calculated free energies 
explicitly consider the entropy of ions and solvent molecules. 
The method is thoroughly described in Refs. 167, 69]; our im- 
plementation and additional terms that must be included in the 
free energy to account for dissociation of weak acid groups are 
described in Appendix B. We neglect spatial variations of seg- 
ment densities in the directions lateral to the surface (and thus 
neglect the effect of ion-ion correlations), and determine the 
variation of densities in the direction normal to the surface. 

All calculations consider the conditions used for the exper- 
imental data shown in Fig. [T^, with pH=5, and 100 mM 1:1 
salt with an additional 5mM 2: 1 salt. The total surface den- 
sity of functionalized molecules (neutral and charged) is ~ 3 



nm ^ at the nanoparticle surface 170(1 . Rather than explicitly 



modeling two species of surface molecules (neutral and acid- 
terminated), we vary the surface charge density ctc by chang- 
ing the valency of the ionic end group. To explore a wide 
range of surface charge we consider end groups with valen- 
cies V G [0,2]; the experimental data in Fig. [T] corresponds 
to the range v G [0, 1]. For v = 1 the number of function- 
alized charge groups corresponds to ~ 80% of the charge on 
a BMV capsid. The lattice size and statistical segment length 



7311 are set to 



for polyelectrolyte and TEG molecules ItiIIt^, 
Z = 0.35 nm (which roughly corresponds to the size of a single 
amino acid and the correlation length of water), and there are 
0.5 charges per polyelectrolyte segment, which corresponds 
to the linear charge density predicted by counterion condensa- 



tion 174 



75 



76, 



7711 . Except for charge, all species are treated 



the spatial distributions of polyelectrolyte, TEG, ions, and wa- 



as chemically identical, with the Flory interaction parameter 
X = and dielectric permittivity e = 80). Free energies 
are insensitive to the introduction of unfavorable interactions 
(X > 0) between the alkane segments and hydrophilic species 
or species-dependent dielectric constants. In the strong cou- 
pling limit, the equilibrium density of positive charge due to 
adsorbed subunits is not sensitive to statistical segment length, 
lattice size, polymer length or charge per segments. The disso- 
ciation constant (pK^) for carboxylate groups is not known in 
the vicinity of the surface functionalized gold particle; we use 
pK^ — 4.5, so that approximately 25% of the acid groups are 
dissociated forpH = 5 at full coverage (see below), consistent 
with Ref. 17811 . We do not explicity consider complexation of 
acid groups, a possibility suggested in that reference. We will 
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FIG. 4: The free energy per area at the capsid surface for the adsorp- 
tion of a concave polymer brush onto a spherical brush with opposite 
charge. For varying surface densities of polyelectrolyte (normalized 
by the density in a complete capsid), psurf, the free energy relative to 
Psurf = is shown for nanoparticle functionalization with weak acid 
groups (pKi — 4.5) at ctc = 3 acid groups/nm^ (■), ac = 1 (A) 
and strong acid groups at ac = 3 (-I-), ac ~ 1 (o). There is 100 mM 
1:1 and 5 mM 2:1 salt andpH = 5. 
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FIG. 5: The equilibrium surface density of adsorbed polyelectrolytes, 
Psurf normalized by the surface density of polyelectrolyte on a core 
with a complete capsid. Predictions are shown for surface function- 
alization with weak acid (■) and strong acid (-I-) groups. The bulk 
subunit concentration is 10 ^M, with other parameters as in Fig.|4] 



also consider experiments in which the terminal group is a 
strong acid, with pK^ = 0.0. 

Free energies and equilibrium concentrations of adsorbed 
polyelectrolytes. The equilibrium surface density of adsorbed 
polyelectrolyte can be calculated by considering a system for 
which there is free exchange between surface adsorbed poly- 
electrolytes and a bath containing polyelectrolyte, salt, and 
solvent. The kinetic theory developed in Section lTVl however. 



requires the free energy for non-equilibrium adsorption den- 
sities. We therefore consider a restricted number of grafted 
molecules (which cannot exchange with the bath), as a func- 
tion of surface charge and density of grafted molecules. For 
each grafting density psurf, the spatial distributions of poly- 
electrolyte, TEG, ions, and solvent are determined by mini- 
mizing the total free energy. The free energy, /sf, is calculated 
from Eqs. lBlll - IB17l Note that the density of polyelectrolytes 
(Fig-EJl segments is nonmonotonic with distance from the cap- 
sid, consistent with studies of capsid assembly around nucleic 
acids I36. 651. 



Free energies for two functionalized charge densities with 
weak and strong acid groups are shown in Fig. 2] and the 
equilibrium adsorption densities, Ps„j.f, which correspond to 
the minimum excess free energy, are shown in Fig. |5] The 
adsorption densities (psurf, Psurf) '^^^ normalized with respect 
to the density of dimer subunits in a complete capsid (0.09 
subunits/nm'^ on the interior surface of the capsid with inner 
radius 8.9 nm). The predicted adsorption densities for weak 
acids (pKn = 4.5) are significantly lower than for the case of 
strong acids because the high local density of negative charges 
shifts the carboxylate dissociation equilibrium; for a carboxy- 
late surface density ac = 3 /nrn^ the average dissociated frac- 
tion shifts from 75% for an isolated molecule to 35%. How- 
ever, dissociation is enhanced as polyelectrolytes with positive 
charges adsorb and the average dissociation fraction on a core 
with a complete capsid is 78%. As discussed in Section [Vl 
this effect results in significantly different dependencies of as- 
sembly kinetics on p^'^^^ for weak and strong acids. 

There are several reasons why our calculations may un- 
derestimate p^Jjj^ at low ac and for weak acid groups. 
Nanoparticle-protein interactions in addition to those involv- 
ing the N-terminal arms could contribute to subunit adsorp- 
tion. The calculated interactions involving N-terminal arms 
assume uniform charge densities in directions parallel to the 
surface, which will overestimate polyelectrolyte interactions 
at low surface densities, when the average distance between 
adsorbed subunits is greater than the polyelectrolyte radius of 
gyration; the approach described for neutral subunits in Ap- 
pendix A might be appropriate in this regime. As discussed 
above, reducing restrictions on polyelectrolyte conformations 
due to the folded portion of capsid proteins results in some- 
what higher adsorption densities for low surface charges. In 
addition, the pK-^ of carboxylate groups and the effect of in- 
teractions between the charged species and the gold surface 



6 




E 



Oq [acid groups / nm 




FIG. 6: The electrostatic potential at the nanoparticle surface for 
varying surface functionalization (crc) with strong (+) and weak (■) 
acid groups. There is 100 mM 1:1 and 5 mM 2:1 salt. 



are poorly understood. For this reason, we examine assembly 
over a wide range of surface charges for both pK^ = 4.5 and 
pK^ = to understand the interplay between surface charge, 
dissociation, and assembly. It would be desirable to experi- 
mentally measure the surface charge as a function of the frac- 
tion of charged TEG molecules on the surface; while it will 
be difficult to quantitatively verify the surface charge or cor- 
responding magnitude of the electrostatic potential, the vari- 
ation of the potential at the core surface with ctc (see Fig.|6]l 
is roughly consistent with preliminary experimental measure- 
ments of zeta potentials for functionalized nanoparticles 
While our calculations may underestimate p^'^^^, these approxi- 
mations are less important at high coverages and hence should 
not significantly affect equilibrium packaging efficiency cal- 
culations. 



A. Empty-capsid free energies. 

We use the same approach, but without the core and TEG 
molecules, to estimate the electrostatic contribution to the free 
energy of forming an empty capsid. The free energy, /sf.ec, as 
a function of polyelectrolyte density (psurf) is shown in Fig.|7] 
where it is compared to the free energy calculated with the 
linearized Poisson Boltzmann equation for charge smeared 
on a spherical surface 1351 16111 with radius i?cap = 8.9 nm, 
/pB = ABADgs^n2/[2i?cap(AD + -Rcap)], with ris = 90 subunits, 
(7s=18 charges/subunit, Ab ~ 0.7 nm the Bjerrum length of 
water, and Ad the Debye length. Although the linearized Pois- 



FIG. 7: The free energy per area at the capsid surface of an empty 
capsid due to electrostatic repulsions and ion and solvent entropy 
is shown as a function of the charge density on capsid subunits, as 
predicted by SF calculations (solid lines with symbols) and the lin- 
earized PB equation (dashed line). The calculations use 100 mM 1:1 
salt and 5 mM 2: 1 salt. 



son Boltzmann equation was shown to agree closely with the 
full nonlinear calculation in Ref. 16111 for this ionic strength, 
there is a large discrepancy between the PB and SF calcu- 
lations; the SF calculations predict much lower free ener- 
gies because the flexible tails adopt stretched conformations 
which reduce electrostatic repulsions (see also Ref. 16511 ). We 
note, though, that we have not considered interactions be- 
tween charges located in the structured region of capsid pro- 
teins. 



Zlotnick and coworkers I29ll33ll79ll have shown that capsid 



formation free energies can be fit using the law of mass action 
to experimental data for capsid and free subunit concentra- 
tions as the total capsid protein concentration is varied. The 
resul ting free energies depend on pH as well as salt concentra- 
tion 113 311 and likely include effects ari sing from hydrophobic 
as well as electrostatic interactions ||35|1 . We can subtract 
the polyelectrolyte free energy from the total capsid forma- 
tion free energy derived from experimental data to estimate 
the free energy due to attractive subunit-subunit interactions, 
gn, which may arise from hydrophobic contacts. Following 
Ceres and coworkers 113311 . we assume that each subunit in a 
capsid makes favorable contacts with n*^ = 4 neighboring sub- 
units, so the attractive energy per subunit is 



5h = <?b«72 - 47ri?2 /sF,Ec/iV 



(1) 



For iV = 90 protein dimer subunits in a T=3 BMV cap- 
sid with an inner capsid radius of i?cap = 8.9 nm 1I57|1 . we 
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obtain /sf.ec = 2.35 ksT/nm^. With a typical subunit- 
subunit contact free energy of = —5 ksT (-3 kcal/mol) 
129, 33 , 79I1 . this gives ~ —36 fcBT/subunit (compare to 



79(1 , this gives ga : 
the estimate from experimental data for hepatitis B virus as- 
sembly by Kegel and van der Schoot ll35ll of an attractive in- 



teraction strength of —13ksT per subunit). 



III. EQUILIBRIUM THEORY FOR CORE-CONTROLLED 
ASSEMBLY 



In this section we consider the thermodynamics for a sys- 
tem of subunits that can assemble into capsids around cores 
and assemble in bulk solution to form empty capsids. The 
free energies due to electrostatic interactions for subunits ad- 
sorbed on surfaces or assembled into empty capsids developed 
in section|II]are used to determine the relative free energies of 
core-associated capsid intermediates. 



A. The thermodynamics of core-controlled assembly 

We consider a dilute solution of capsid protein subunits 
with density pp, and solid cores with density pc- Subunits can 
associate to form capsid intermediates in bulk solution or on 
core surfaces. A complete capsid is comprised of N subunits. 
Zandi and van der Schoot [? ] develop the free energy for 
core-controlled assembly in the context of capsids assembling 
around polymers, but neglect the contribution of partial cap- 
sid intermediates, which will be important for describing the 
kinetics of core-controlled assembly in the next section. Here, 
we consider the free energy for a system of subunits, capsids, 
and partial capsid intermediates. To simplify the calculation, 
we neglect the possibility of malformed capsids, but note that 
disordered configurations could be kinetically arrested or even 
the favored equilibrium state in the case of core-subunit inter- 
actions that are much stronger than the thermal energy k^T 
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The total free energy density is given by 



F = Fec + Fa 



(2) 



with Fec the free energy density of empty capsid intermedi- 
ates 



N 



EC 



Pi log{pia ) - Pi+ PiG.^ 



cap 



(3) 



where pi is the density of intermediates with i subunits and 
a? is a standard state volume. For simplicity we follow Zlot- 
nick and coworkers ll3(Xl3lll and assume that there is only one 
intermediate for each size i, with a free energy Gf^ defined 
below in Eq. |9] 

The free energy density of core-associated intermediates is 
given by 



F„ 



PC 



N N 



pea ) 



n— m— 



n,m r -i n,m^n,m 



(4) 



where Pn,m is the fraction of cores with an adsorbed inter- 
mediate of size m and with n adsorbed but unassembled sub- 
units. The free energy for such a core is given by G™^,jj and 
is defined below. For simplicity, we assume that a core can 
have at most one assembled intermediate and a total of N sub- 
units can fit on a single core surface, where N is the size of 
a complete capsid, so states with m + n > N subunits have 
zero probability. The prime on the second sum indicates that 
TO = 0, 2, 3 ... because, by definition, a state with m = 1 
has no assembled intermediates. 

To obtain the equilibrium concentration of intermediates we 
minimize the total free energy subject to the constraints 



N N ' 
n— 1 rn— 



N 

Ew 



Pv 



(5) 



and 



N N 



E E ~ ^ 



(6) 



where Eq.|5]enforces that the total density of subunits is given 

by Pp- 

Minimization yields 



(7) 



and 



i=l 



Pn,mPca^ = exp[-/3(G^„°;^ - {n + m)p - ^c)] (8) 

where j3 = I/Zcb^ is the inverse of the thermal energy, 
p ~ kBT\og{a^ pi) and pc = kBT\og{a^ pcP^fi) are chemi- 
cal potentials for subunits and cores, respectively. 
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B. Free energies for empty-capsid assembly intermediates 

The free energy of an empty capsid intermediate with i sub- 
units is written as 



ksTS', 



fdegen 



iksTsrot (9) 



where tij is the number of new subunit-subunit contacts 
formed by the binding of subunit j to the intermediate, g\, 
is the contact free energy, S'^'^^™ accounts for degeneracy in 
the number of ways subunits can bind to or unbind from an 
intermediate (see the s factors in Refs. IsO, 31 1), and Srot is 
a rotational binding entropy penalty. For simplicity and to 
ease comparison with earlier works, we take Srot = except 
when comparing theoretical predictions to simulation results 
(see Appendix A). As discussed at the end of Section |ll] the 
subunit-subunit contact free energy has been estimated from 



experiments for several viruses [29 
on salt concentration and pH 133 



33 



59 



3511 



7911 and depends 



C. Free energies for assembly intermediates on core surfaces 

The free energy for a core without a partial capsid is given 
by G™'q= = Ac f SF {(Jc, P>iuvf) with Ac the core area, ctc the 
functionalized acid group density, psurf = n/N the normal- 
ized density of adsorbed subunits, and /sf calculated in Sec- 
tion ini Likewise, the free energy for a complete capsid is 
6™^"= = G'J + Ac[/sF(ac, 1) - /sF,Ec], where the polyelec- 
trolyte contribution to the capsid free energy (/sf.ec) is sub- 
tracted because it is included in /sf- 

In the case of a core with a partially assembled capsid and 
adsorbed but unassembled subunits, there is an inhomoge- 
neous distribution of charge on the surface, with higher charge 
density in the region of the assembled cluster We determine 
the free energy for such a core, with n adsorbed subunits and 
an intermediate of size m, in the following schematic way. 
We assume the unassembled subunits spread on the surface 
of the core not occupied by the intermediate. There are thus 
two 'pools' of charge on the surface - the partial capsid with 
a high subunit density psurf ~ 1 and the remaining area of the 
core surface with a charge density spread uniformly over the 
area not covered by the intermediate , A^*^ = Ac {N — m)/N. 
The total electrostatic free energy is then given by 



GT,L = + {Ac - A^^T) [/sF(ac, 1) - /SF,EC] 



^mVsF(o-C, 



00 



nqp 



N 



- ) for n + m < N 

for n + ni>N. (10) 



In general there should be a term in G"^™'' to account for strain 
energy if the core curvature is not commensurate with that of 
the intermediate; we assume that capsid and core curvatures 
match. 

Core encapsidation. Above a threshold or critical sub- 
unit concentration (CSC), the majority of subunits will be in 
capsids[3C , 
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SOU . Because core-subunit interactions can 



lead to high localized surface charge densities, the thresh- 
old concentration for assembly in the nanoparticle system 
can be below the CSC for spontaneous assembly. Above 
the CSC, both empty and core-filled capsids can assemble. 
Due to the entropy cost of assembling a capsid on a core, 
Eqs.[8]and|2]show that there is a threshold surface free energy 
(G™'^Q — G^'' « — feB^ log(pca'^) for a stoichiometric amount 
of capsid protein and nanoparticles), below which empty cap- 
sids and free cores are favored over full capsids. Above the 
threshold surface free energy and above the CSC, nearly all 
cores will be encapsidated (Po.iv ~ 1) at equiUbrium. How- 
ever, we will see that there can be a higher threshold surface 
free energy for efficient encapsidation kinetics. 



IV. KINETIC THEORY FOR CORE-CONTROLLED 
ASSEMBLY 



In this section we use the free energies developed in sec- 
tion|lII]for capsids and capsid intermediates as the basis for a 
kinetic theory of core-controlled assembly. 

Zlotnick and coworkers BOl 13111 have developed a system 
of rate equations that describe the time evolution of concen- 
trations of empty capsid intermediates 



dpi 
dt 



— = -2fcasS2Pl 



N 



dpi 
~dt 



kzsSiPlPi-l — fcas^i+lPlPi 
-k^Pi + fcbPi+l 



.N 
(11) 



where fc^s and k^ are respectively the forward and reverse 
binding rates and is a statistical factor that describes the 
number of different ways to transition from intermediate i to 
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i + 1 i3(X blh ; S2 corrects for double counting the number of 
pairwise combinations of free subunits. Following Ref. h31I1 . 
we simplify the calculation by requiring that transitions be- 
tween intermediates are only allowed through binding or un- 
binding of a single subunit and there is only one intermediate 
for each size i. We assume that the assembly rate constant 
fcas is independent of intermediate size i, but the number of 
contacts per subunit is not (see Table 1). 

We extend this approach to describe simultaneous capsid 
assembly and adsorption to spherical cores, giving the follow- 
ing expression for the time evolution of core states 



dt 



-Wr^-2.2Pn-2,2+WT"'Pna 



dPn: 
dt 



^n-l,2PlPn-1.2 + 2^71+1, 2 



— {plVLn,2 + ^n,2^ Pn,2 

+ Wn+2fiP7i+2fl + W^n-l,3^n-l,3 
+ Wr'Pn,0+WT"Pn,3 

- (Wn,2 + W„,2 + PlM/r" + WT') Pn,2 
dPn. 



dt 



— Pl^n—l,7n^n—l, 



m^3...N 



,m+l Pn-l ,m+l 
+PlW^_^Pn^ni-l + iy„+iP„,,„+l 

- f W„,™ + W„.™ + Pl Wr"' + "^^r"') Pn.n (12) 



The first two lines in each of Eqs. [12] describe adsorption 
and desorption from the core surface with respective rate con- 
stants and n, while the following lines describes binding 
and unbinding to a capsid intermediate. Rate constants for 
binding and unbinding of adsorbed subunits to an adsorbed 
intermediate are denoted by W and W, respectively. Finally, 
]/l/spont ^jj^ yy*"''""' are respectively the rate constants for the 
association (dissociation) of non-adsorbed subunits to (from) 
an adsorbed intermediate. The latter process becomes im- 
portant when an adsorbed capsid nears completion, and elec- 
trostatic repulsions render non-specific adsorption of subunits 
unfavorable, causing the surface capture and diffusion mech- 
anism to be slow in comparison to the empty-capsid assembly 
mechanism. 



Adsorption is calculated from Q.n,m = fcad'&n.m, with fcad 
the adsorption rate constant and the probability that an 
adsorbing subunit is not blocked by previously adsorbed sub- 
units. For simplicity, we assume Langmuir kinetics, ^n,m = 
{N — m — n)/N, but note that this is a poor approximation at 
high surface coverages (see Appendix A). 

The desorbtion rate ri„ m is related to the adsorption rate 
by detailed balance 

a„„ = r!„_i.™ exp {Gl"% - GT-\J /a\ (13) 

Assembly and disassembly of adsorbed subunits is de- 
scribed in a manner analogous to that used for empty capsids, 
with assembly rates given by 



(14) 



surf 



Dr 



m^2...N-l 



where the effective surface density is p^^^ — n/[a^{N — m)]. 

Eq. [15] states that subunit-subunit binding rates are propor- 
tional to the frequency of subunit-subunit collisions, which 
can be calculated from the Smoluchowski equation for the 
diffusion limited rate in three dimensions, with the density of 
subunits iPn'^) within a layer above the surface with thick- 
ness of the subunit size, a, or from the diffusion limited 
rate in two dimensions 18 ill with a surface density given by 
P2D — Pn^m^- The rcsult of the two calculations differs only 
by a logarithmic factor that is of order 1 in this case (see 
appendix B of Ref. Isil]). Surface assembly rate constants 
(/cas) can be smaller than those for empty capsid assembly 
(/cas), if interactions between subunits and the TEG layer im- 
pede diffusion or reorientation of adsorbed subunits. Associ- 
ation rates could also change in comparison to bulk solution 
if subunits take different conformations on the core surface; 
however, imaging experiments show that protein structures 
in nanoparticle-filled capsids are nearly identical to those in 
empty capsids i24ll . which suggests that subunits do not dena- 
ture on the surface. 

The disassembly rate is given by detailed balance 



Wr. 



W^+,,m-i exp (GTm - /a\ (15) 



Similarly the association and dissociation rates for solubi- 
lized subunits are related by detailed balance: 



^spoilt 



hptsS 



.N - 1 



w: 



■spont 



W^r" exp {GZl ~ _i) /a' (16) 
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The expression for the time evolution of empty capsid in- 
termediates is now 



(a) 1 



dpi 
dt 



N 



-2fcasS2/Oi - ^ hsSi+lPiPl + hpt 



i=2 



N N 



PC 



n^,m+W'r'-piWr')Pn,. 

i = 2..N 



dp. 



-TT = KsStPlPi-1 - KsSi+lPlP, 

at 

-k^Ui + k^Pi+i. 



(17) 



We note that although Eqs. [12] and [17] have the form of a 
Master equation, the finite concentration of subunits intro- 
duces a nonlinear dependence of the time evolution on state 
probabilities. 

V. RESULTS 

In this section we report the predictions of the equilib- 
rium and kinetic theories for two experimentally observable 
quantities: packaging efficiencies, or the fraction of nanopar- 
ticles encapsulated by well-formed capsids (estimated from 
TEM experiments), and time-resolved light scattering inten- 
sity. Preliminary experimental data for these quantities is 
shown in Fig. [T] Since there is not sufficient data to esti- 
mate the values of all unknown parameters, we make qualita- 
tive predictions about the effect of changing the experimental 
control parameters: the functionalized surface charge density 
CTc, the stoichiometric ratio of capsid protein to nanoparticles 
''CPiNP, and the capsid protein concentration Cc- The impor- 
tant unknown parameters are the subunit-subunit binding en- 
ergy, gb, the dissociation constant of the functionalized charge 
groups, the rate constant for adsorption of subunits onto the 
nanoparticle surface, fcad, and the rate constant for the assem- 
bly of surface-adsorbed subunits, fcas- The definitions and val- 
ues of model parameters are given in Table |T] 

A. Packaging efficiencies 

We begin by discussing the relationship between the density 
of functionahzed charge groups on nanoparticle surfaces, ctc, 
and packaging efficiencies. As discussed in Section|l] the ex- 
perimental measurements of packaging efficiencies were ob- 
tained from an experimental protocol in which the solution pH 
varied over time; because we do not know the dependence of 
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FIG. 8: Packaging efficiencies, or the fraction of nanoparticles incor- 
porated in well-formed capsids, at varying surface charge densities 
((Tc) of weak acid groups predicted by (a) the equilibrium theory 
and (b) the kinetic theory at 10'' seconds. Parameters are as listed 
in Fig.Ewith Ka = 10^(M ■ s)"\ = 10*M s"\ and subunit 
binding free energies are indicated on the plots. 



the subunit-subunit binding energy on pH, we only present 
predictions of the equilibrium and kinetic theories for packag- 
ing efficiencies at pH = 5 for several binding energies. We 
will compare these predictions to the currently existing data, 
but note that they should be compared to packaging efficien- 
cies from experiments carried out entirely at pW = 5, as was 
the case for the light scattering data shown in Fig.[T}3. 

As shown in Fig. [8^, the equilibrium theory predicts effec- 
tive packaging only above a threshold functionalized charge 
density, which depends on the binding free energy. Effec- 
tive packaging even occurs at a binding free energy — 0, 
showing that assembly on cores can occur under conditions 
for which spontaneous empty-capsid assembly is not favor- 
able, as seen in experiments [24.] and simulations ll43ll44|] . On 
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core surfaces nanoparticle-subunit electrostatic interactions 
and attractive subunit-subunit interactions overcome subunit- 
subunit electrostatic repulsions to drive assembly (gu = 
— 15.6 kcal/mol per subunit for g\, ~ 0, Eq.[T]i. 

As shown in Fig. [8j), for weak subunit-subunit interactions 
(gb > —2 kcal/mol) the kinetic theory predictions at 10^ sec- 
onds are essentially identical to the equilibrium results. As 
subunit-subunit interactions increase in magnitude, however, 
the predicted threshold charge density becomes larger than the 
equilibrium value, and even begins to increase for < —3. 
With stronger subunit-subunit binding interactions, the assem- 
bly of empty capsids competes with assembly on core sur- 
faces by depleting the concentration of free subunits. Assem- 
bly on cores remains favorable at high surface charge den- 
sities because initial subunit adsorption leads to rapid nucle- 
ation, and subsequent nonspecific adsorption enhances capsid 
growth rates (see below). The predicted packaging efficien- 
cies at long times are relatively insensitive to variations of the 
kinetic parameters fcad and fcas. 

Interestingly, the threshold surface charge density predicted 
by the kinetic theory, uc ~ 1-2, is nearly independent of the 
subunit binding energy for gb < 2 kcal/mol, and comparable 
to the lowest experimental charge densities with measurable 
packaging efficiencies (Fig. [U. However, the experimental 
results show a gradual increase in packaging efficiency as the 
functionalized charge density increases. We do not observe 
such a gradual increase in packaging efficiencies except for 
shorter observation times or an extremely low surface assem- 
bly rate constant. This discrepancy could arise because we 
have assumed monodisperse cores that are perfectly matched 
to the capsid size, or it may indicate the presence of kinetically 
frustrated states (see below). 

B. Assembly kinetics on core surfaces 
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In this section we explore the effect cargo-capsid protein 
interaction strength on assembly kinetics. As discussed in 
Section U assembly kinetics can be monitored in experiments 
with time-resolved dynamic light scattering; preliminary ex- 
perimental data for the highest functionalized charge density 
is shown in Fig.[T] Because it is unlikely that the light scatter- 
ing signal can distinguish between disordered and assembled 
subunits on core surfaces, we estimate fight scattering from 
the kinetic theory by calculating the mass-averaged amount 
of adsorbed and assembled (n + m) subunits on core surfaces. 



FIG. 9: The predicted time dependence of light scattering for func- 
tionalization with strong acid groups, (a) and (b) Indicated func- 
tionalization densities with rate constants (a) fcas = lO'^M s~^ and 
(b) fcas = lO^Ms"^. (c) Indicated surface assembly rate con- 
stants, with (Tc = 3 nm~^. Other parameters for (a)-(c) are fcad = 
\'d'{M s)"\rcP:NP = 1.33, Cp = lO^iM, and = -2 kcal/mol. 

We simplify the presentation by varying only the charge den- 
sity (ctc) and the surface assembly rate constant (/cas), with 
physically reasonable values for the subunit binding free en- 
ergy, 5b = ~2 kcal/mol, and the the subunit adsorption rate 
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FIG. 10: The predicted time dependence of ligiit scattering for sur- 
face cliarge groups with pK^ — 4.5 and pH = 5, for (a) indicated 
functionalization densities fcas = lO'^M s~^, (b) indicated function- 
alization densities with fcas = lO^M s~^, and (c) indicated surface 
assembly rate constants with ac = 3 nm^^. Other parameters are as 
inFig.m 



constant, fcad ~ 10^(M ■ s)~^ ll82ll . As discussed in section 
Section IIVI the effect of reduction in dimensionality due to 
adsorption onto a surface is accounted for in the kinetic the- 
ory, but assembly rate constants could still be dramatically 



different from those in bulk because of impeded diffusion of 
adsorbed subunits. We therefore explore a wide range of sur- 
face assembly rate constants fcas- 

Predicted light scatter as a function of time is shown for 
various functionalized charge densities ctc and surface assem- 
bly rate constants fc^s in Figs^ and [TO] for functionalization 
with strong and weak acids BSSll . respectively. The analysis is 
easiest for the case presented in Fig.|9^ (fcas = lO^M s~^), for 
which assembly on the core surface is slow compared to sub- 
unit adsorption. For each value of ac, there is a rapid initial 
rise in light scatter, followed by a charge-dependent plateau, 
and finally a slow increase to saturation. The initial fast ki- 
netics correspond to nonspecific subunit adsorption (i.e. with 
little or no assembly) and the plateau occurs when adsorption 
saturates at what would be the equilibrium surface density if 
there were no assembly ipl'^^, see Section|ll|l. If adsorbed sub- 
units assemble into a relatively stable intermediate, the chem- 
ical potential of remaining adsorbed subunits decreases, re- 
sulting in further adsorption and increased light scatter. This 
process continues as additional subunits bind to the growing 
intermediate, until assembly stops. 

At high functionalization densities, nonspecific subunit ad- 
sorption approaches the density of a complete capsid; the high 
local or surface concentration of subunits enables rapid nu- 
cleation and assembly and thus the light scatter reaches its 
final saturation value at relatively short times. At the lower 
functionalization densities, nonspecific adsorption saturates 
at smaller surface concentrations, with corresponding smaller 
initial plateau heights for the light scatter. Nucleation times 
increase dramatically as the surface concentration is reduced 
(see the discussion in Ref. 14411 '). resulting in the long plateau 
before assembly completes. 

For a larger assembly rate constant fcas = 10''M s^^, 
the nonspecific adsorption and assembly phases can only be 
clearly distinguished at low functionalized charge densities, 
as shown in Fig. |9j3; at larger ac adsorption and assembly 
are concomitant. For weak acid groups (Fig.fToli. equilibrium 
surface densities increase only slowly with ac ipl'^^ < 0.15, 
see Fig.|5]) because the local concentration of negative charge 
disfavors carboxylate dissociation. There is still a strong de- 
pendence of overall assembly rates on ac, however, because 
the positive charge from capsid proteins drives additional car- 
boxylate dissociation (see Section Thus, assembly effec- 
tiveness cannot be predicted from pl'^^^ alone, as was found for 



simulations of neutral subunits in Ref. ll43ll . 
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Cores can enhance assembly rates. From the preceding 
discussion, it is clear that functionalized charge leads to an 
increased concentration of capsid proteins nanoparticle sur- 
faces relative to that in bulk solution, which can increase 
capsid nucleation rates and even induce assembly below the 
threshold protein concentration for formation of empty cap- 
sids. Because nonspecific adsorption onto core surfaces can 
be much faster than subunit-subunit binding rates, the core 
surface can also enhance assembly rates during the growth 
phase by increasing the effective subunit-subunit interaction 
cross-section. This effect was hypothesized for RNA by Hu et 
al. II42II and is reported for simulation results in Ref. ll4411 . 

Assembly rates decrease as capsids near completion. Net 
assembly rates decrease as capsids near completion because 
adsorption of subunits onto cores with large partial capsids 
is impeded by two effects. First, the excluded volume of the 
partial capsid decreases adsorption rates; a similar effect has 
been observed in simulations representing empty capsid as- 
sembly, but the magnitude of rate reduction seems to be model 
dependent (see Refs. 144)^ 149!] ). Electrostatic repulsions with 
adsorbed and/or assembled subunits on the core surface can 
have a larger effect on subunit adsorption. When the subunit 
surface density is large compared to p^^^, the surface capture 
and diffusion assembly mechanism that enables rapid core- 
controlled assembly becomes inefficient. Even though direct 
association of free subunits onto core-bound intermediates is 
included in the kinetic theory, there is a reduction in assembly 
rates as capsids near completion; this effect is most noticeable 
for surface functionalization with weak acid groups (Fig.fTok). 
or low fractions of strong acid groups (tJc < 2 nm^^). The 
decreased assembly rates are reflected in slow growth of light 
scatter as assembly nears completion. 

Comparison with experimental data. As discussed above, 
there is not sufficient experimental data to estimate values for 
the unknown parameters; however, the predicted light scatter 
plots share some features with the experimental light scatter 
data shown in Fig. [TJ). In particular, for some parameter sets 
light scatter increases rapidly at short tirnes, without the lag 
time seen for empty capsid assembly ifssllsg . 601. It is tempt- 
ing to attribute the initial rise to nonspecific adsorption that 
does not saturate until nearly complete coverage, as seen for 
strong acid surface groups (Fig. |9^), but the calculation with 
weak acid groups predicts values of p^^^ that are well below 
complete coverage even for ac = i nm~^. Interestingly, the 
predicted light scatter curves in this case with fast assembly 
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FIG. 1 1 : Estimated light scatter from simulations in which subunits 
assemble around model nanoparticle with a size mismatch. The 
preferred empty-capsid morphology is a T=3 capsid, but the low- 
est free energy configuration is a T=l capsid assembled around the 
T=l-size nanoparticle. Estimated light scatter is shown as a function 
of time for different values of the parameter eaa (in units of fcBT), 
which gives the free energy cost for adopting a subunit conformation 
consistent with a T=l morphology. As described in Ref. l43ll , the 
larger values of eaa lead to the formation of frustrated nanoparticle- 
associated partial capsids with a T=3 morphology, which cannot 
close around the nanoparticle and block assembly into the lowest free 
energy configuration. The data is replotted from simulations reported 
in Fig. 6 of Ref. l43ll . with the nanoparticle-subunit interaction en- 
ergy = S/cBr. 



kinetics (fcas — lO^M s^^, Fig.|9^) still demonstrate a rapid 
initial rise followed by slow saturation; the initial rise coiTe- 
sponds to rapid nucleation and assembly concomitant with ad- 
sorption. Thus, as discussed in the next section, comparison 
of theoretical predictions and experimental light scatter data 
at varying functionalized charge density will be necessary to 
unequivocaUy determine the origin of the initial light scatter 
phase. 

The experimental light scattering data (Fig. [T}^) appears to 
demonstrate a logarithmic increase in light scattering between 
20 and 200 seconds. While the theoretical light scatter curves 
also demonstrate slow growth near completion (see above) 
due to subunit-subunit electrostatic repulsions, this effect is 
limited to the last five or 10 subunits. We did not find any 
parameter sets for which predicted light scatter grew logarith- 
mically over several decades in time, although our parameter 
search was not exhaustive. Intriguingly, the predicted light 
scatter from simulations in which frustrated off-pathway in- 
termediates block assembly ll43ll do show logarithmic growth 
over several decades (see Fig.fTTTl. 
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C. Concentration effects 

Varying nanoparticle concentrations. We first examine 
assembly effectiveness as the nanoparticle concentration is 
varied at fixed subunit concentrations. We define the cap- 
sid protein-nanoparticle stoichiometric ratio as rcp:NP = 
Cp/{NCc) so that rcp:NP = 1 vvhen there is exactly enough 
capsid protein to assemble a complete capsid on every core. 
The variation of packaging efficiencies with the stoichiomet- 
ric ratio is shown in Fig. [121 the predicted curves are sigmoidal 
in shape, much like the experimentally observed packaging ef- 
ficiencies shown in Fig. lA of Ref. l24ll . Packaging efficien- 
cies for rcp:NP ^ 1 are typically lower than the equilibrium 
values obtained by solving Eqs.|2]and[8] because initially cap- 
sid proteins undergo nonspecific adsorption onto every core, 
which depletes the concentration of free subunits. The forma- 
tion of complete capsids then requires desorption and subse- 
quent re-adsorption onto cores with large intermediates. The 
timescale for subunit desorption increases with functionalized 
surface charge density and therefore the kinetic packaging ef- 
ficiency is nonmonotonic with respect to ac for rcp:NP < 1- 
This scenario is consistent with the slow assembly kinetics ex- 
perimentally observed for capsid assembly in the presence of 
RNAIUJ]. 

Varying subunit concentrations. As shown in Fig. 1131 at 
(Tc > 2 acid groups/nm^ the kinetic theory predicts a much 
weaker dependence of assembly effectiveness on subunit con- 
centration than has been seen for empty-capsid assembly. This 
trend arises because the chemical potential of adsorbed sub- 
units is dominated by electrostatics rather than subunit trans- 
lational entropy. 

VI. DISCUSSION 

The feasibility of estimating the unknown parameters. 

The results in the last section are consistent with some features 
of the preliminary experimental data, and the predicted de- 
pendencies of packaging efficiency and light scattering on the 
functionalized charge density and nanoparticle-protein stoi- 
chiometric ratio can be qualitatively compared with additional 
experimental data. Quantitative comparisons of experimental 
data, however, will require estimates for at least three cur- 
rently unknown parameters, the subunit-subunit binding en- 
ergy g\,, the assembly rate constant on core surfaces fc^s, the 
subunit adsorption rate k.^^. In Appendix A, we show that. 
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FIG. 12: Packaging efficiencies depend on the capsid protein/core 
stoichiometric ratio. The predictions of the kinetic theory for the 
variation of packaging efficiencies with the capsid protein/core sto- 
ichiometric ratio, rcp:NP = Cp/NCc at a fixed subunit concentra- 
tion of Cp — 10/iM, are shown for functionalization with weak acid 
groups, with other parameters as in Fig. [8] 
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FIG. 13: Packaging efficiencies depend only weakly on the initial 
subunit concentration, Cp, for fixed rcpiNP. The predictions of the 
kinetic theory for the variation of packaging efficiencies with sub- 
unit concentration Cp are shown for function, 
groups, with other parameters as in Fig. 1 121 



with good estimates for the unknown parameters, the theoret- 
ical predictions agree reasonably with simulation results for a 
model of core-controlled assembly. While it is not certain that 
a simplified theory with a few parameters will exactly match 
experimental data from a system with thousands of distinct in- 
termediates, each of which could have different rate constants, 
rate equations with a similar level of coarse graining have suc- 
cessfully explained many features of empty-capsid assembly 
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Given the number of experi- 
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mental control parameters in the nanoparticle-capsid assem- 
bly system, it should be possible to generate sufficient experi- 
mental data to further test the qualitative trends predicted here, 
and to make quantitative estimates for the unknown parame- 
ters through data fitting. With estimated values for fcas and fc^d, 
it would be possible to predict the assembly state of subunits 
on core surfaces during the initial rise in light scatter, which 
would be challenging to determine with experiments alone. 

In addition to measuring packaging efficiencies and light 
scattering data over a range of functionalized charge densities 
(with both weak and strong acid groups) and capsid protein- 
nanoparticle stoichiometric ratios, the theoretical predictions 
in the last section suggest additional experiments that could 
be useful to provide independent estimates of some unknown 
parameters. In particular, the predicted light scatter curves in 
Figs. |9] and [TOl show that adsorption and assembly happen si- 
multaneously for many parameter sets. Measuring light scat- 
ter of assembly incompetent proteins in the presence of func- 
tionalized nanoparticles would enable independent estimation 
of the subunit adsorption rate and the equilibrium subunit sur- 
face charge density p^'^^^ ( Fig. |5]). In addition to varying the 
functionalized charge density, the subunit-nanoparticle inter- 
action could modified by mutations to the N-terminal arm on 
capsid proteins ll84ll : Tang et al. 18511 have shown that cowpea 
chlorotic mottle virus capsid proteins with 34 residues deleted 
from the N-terminal arms are assembly competent (although 
they lose selectivity for the native capsid structure). 

Identifying metastable disordered states. We make sev- 
eral important approximations in this work. Most signifi- 
cantly, the theory considers only only intermediates consis- 
tent with the native capsid geometry, while experiments BTOll 
and simulations 143115511 indicate that other capsid morpholo- 
gies and/or asymmetric, malformed structures can be kineti- 
cally or even thermodynamically favored when core-subunit 
interaction strengths are large compared to the thermal en- 
ergy k^T. These effects could be accounted for in the theory 
by extending the state space to include structures other than 
complete and partial well-formed capsids, and by introduc- 
ing subunit diffusion constants that depend on the strength of 
core-subunit interactions. If the present theory proves capable 
of describing experiments with parameters for which predom- 
inately well-formed capsids assemble, however, inconsisten- 
cies between theoretical predictions and experimental mea- 
surements for other system parameter values could identify 
kinetic traps. For example, the calculated packaging efficien- 



cies show a less gradual increase with functionalized charge 
density than measured in experiments and predicted light scat- 
ter does not fully capture the a logarithmic increase in the 
experimental light scatter curve. Both of these discrepancies 
could be explained by frustrated off-pathway of assembly in- 
termediates or slow subunit diffusion rates, as evidenced by 
predicted light scatter from simulations which do demonstrate 
frustrated assembly (Fig.fTTI). Confirming this possibility will 
require additional experimental data. 

VII. CONCLUSIONS 

In summary, we develop simplified thermodynamic and ki- 
netic theories that describe the simultaneous assembly of vi- 
ral proteins and encapsidation of cargo. When applied to the 
encapsidation of charge-functionalized nanoparticles, the the- 
ories predicts a transition from inefficient core encapsidation 
to nearly 100% incorporation efficiency as the functionalized 
charge density is increased beyond a threshold value, and the 
estimated light scatter signal shows an initial rapid increase 
followed by slow rise to saturation, as seen in experiments. 
The predicted increase in incorporation efficiency with in- 
creasing charge density is sharper than seen in experiments, 
which could indicate the presence of kinetic traps that are not 
accounted for in the present theory; comparison of theory pre- 
dictions with experimental data collected with a wide range 
of control parameters will be important to assess this possi- 
bility. The trends predicted here for varying surface charge, 
subunit concentration, and capsid-nanoparticle stoichiometric 
ratio could guide the design of experiments that identify fun- 
damental principles and/or additional complexities for simul- 
taneous assembly and cargo encapsulation. 
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APPENDIX A 

In this appendix we compare predictions of the kinetic the- 
ory for the time dependence of the number of adsorbed parti- 
cles, n + m, to results presented in Ref. [44,1 from a computa- 
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TABLE I: Parameter values used for calculations in this work. The parameters used for capsid assembly are as follows. The number of 
subunit-subunit contacts are ri^ — 1 for j < Unuc, rfj — 2 for j £ (n„uc, N — Unuc + 1], Ji-jv-i = 3 for J ^ — nnuc,N), and n% — 4. This 
choice obviates the need for different nucleation and elongation rate constants (Slll . The nucleus size is rinuc = 5. The forward and reverse 
reaction degeneracies are S2 = S2 ~ 2, Sj = Sj = S for 3 < j < N (the average value calculated from simulations in Ref. (3]) and S]v = 1, 
sjv = A^. The binding degeneracy entropy (Eq.[9} is 5"''=s'=" — f^^ log{sj /sj). 



tional model that represents the assembly of T=l capsid-like 
objects around a commensurate nanoparticle. This compari- 
son elucidates the effect of approximations used in the kinetic 
theory. In particular, in rate equation approaches that assume 
a single reaction pathway as described above, or multiple re- 



action pathwa)^ 5C , 5lJ,|87t], and even binding between large 



intermediates 11501 Dill , only intermediates that are consistent 
with partially assembled capsids are considered. The simu- 
lations, on the other hand, explicitly track the spatial coordi- 
nates of subunits and thereby make no a priori assumptions 
about assembly pathways. Because we can either calculate 
parameters or estimate them from independent simulations, 
this comparison is a stringent test of the validity of these ap- 
proximations. It also tests the utility of using the kinetic the- 
ory to describe experimental adsorption results. The proce- 
dure by which we determine parameters for the kinetic theory 
that correspond to particular sets of simulation parameters is 
described below, as well as modifications to the theory in or- 
der to represent neutral subunits (or high salt concentrations). 
We do not consider here the simulations reported in Fig. [TT] 
and Ref. 14311 . for which the nanoparticle is not commensurate 
with the preferred capsid morphology. 

Adsorption kinetics predicted by the kinetic theory and ob- 
served in simulations are shown in Fig. [14] for several subunit 
concentrations and subunit-subunit binding energies for an ad- 
sorption energy of = 7. The agreement between the kinetic 
theory and simulation results is surprisingly good, consider- 
ing that no parameters were adjusted to fit this data and, as 



discussed below, relatively crude estimates are used for the 
subunit-subunit binding rate constants / and binding entropy 

The computational model considered in Ref. consists 
of rigid subunits, for which excluded volume interactions are 
modeled by spherically symmetric repulsive forces, and com- 
plementary subunit-subunit interactions that drive assembly 
are modeled by directional attractions. The lowest energy 
states in the model correspond to "capsids", which consist of 
multiples of 60 monomers in a shell with icosahedral sym- 
metry. The parameters of the model are the energy associ- 
ated with the attractive potential, Sb , which is measured in 
units of the thermal energy /cbT, and the specificity of the 
directional attractions, which is controlled by the angular pa- 
rameters djji and 0,1,. Capsid subunits also experience short 
ranged isotropic interactions with a rigid sphere placed at the 
center of the simulation box; these interactions are minimized 
when capsid subunits are adjacent to the surface of the sphere. 
Subunit positions and orientations are propagated according 
to overdamped Brownian dynamics, with the unit of time 
to = a^/ASD, where D is the subunit diffusion coefficient. 
Full details of the model are given in Ref. 14411 

Free energies for neutral subunits with no internal struc- 
ture or high salt concentrations. The free energies for partial 
capsids on core surfaces are modified from the expressions 
given in Section ITlI CI for the case of neutral subunits or high 
salt concentrations, in which case interactions are limited to 
excluded volume, directional attractions that drive assembly. 
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FIG. 14: The kinetic theory predictions (smooth lines) and simula- 
tion results (noisy lines) for the time dependence of adsorbed and 
assembled subunits, n + m, are shown for the neutral subunit model 
(Eq. lAlb . The simulations consider T=l capsids and commensurate 
cores, so a complete capsid is comprised of Ai' = 30 dimer subunits. 
Curves at increasing height correspond to reduced subunit densities 
of ppo^ = 2.04, 4.07, 8.14, 20.4, 40.7, with a surface free energy of 
Ec = 7 and subunit-subunit binding energies of (a) Sb = 9.0 and (b) 
£b = 11.0 



and favorable core-subunit interactions. The free energy of 
a core with n adsorbed but unassembled subunits and an ad- 
sorbed intermediate of size m is written as 



G: 



surf 



(n + m)e 



TS'. 



mix 

m.n 



TS'"^n, m 



(Al) 



where eh is the core-subunit surface free energy strength. 

The second term, S*™" accounts for the entropy for two- 
dimensional motions of adsorbed subunits on the core surface. 
A simple approach would assume Langmuir adsorption, but 
we find that much more accurate results are obtained by in- 
tegrating an empirical formula for the chemical potential of a 
fluid of hard disks LSSil because the Langmuir model overesti- 
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with the packing fraction 77 = (n + 7n)a^/[167r(i?c + 

The last term in Eq. lAll S'"^^, is the entropy penalty for sub- 
units to be localized a smaller distance from the core surface 
than the subunit diameter, which we obtained by numerically 
integrating the partition function for an adsorbed subunit when 
comparing the kinetic theory to simulations. 

Calculating the time dependence of adsorbed and assem- 
bled subunits requires choosing values for several theoretical 
parameters; here we describe the parameter values and how 
they were chosen. 

1. Subunit binding entropy penalty. Eq.|9] includes a term 
Srot for the entropy loss (in addition to the subunit mix- 
ing entropy included in gb) for a subunit to bind to a cap- 
sid that accounts for translational restrictions on scales 
smaller than the subunit diameter a and rotational re- 
strictions ; the entropy penalty should increase only 
slightly when a subunit changes from one to more than 
one bond. The subunit-subunit binding entropy penalty 
is calculated as described in Ref. 1I45I1 . There is a typo 
in Eq. (15) of that reference; the correct formula is 



3 , /352u,„(r) 



r=0 



1 [Pe^fn^ 
-2^°^^^^ ^^^^ 



2. Subunit binding rate constant /. The assembly rate 
constant was chosen as / = 0.03 (in the dimension- 
less units of Ref. ||44l] ) by comparing kinetic theory 
predictions to simulation results for the assembly of 
empty capsids. The resulting fit (not shown) shows less 
agreement between theory and simulation results than 
we find for core simulations, perhaps because approxi- 
mations such as an average assembly pathway without 
multimer binding are less significant when subunits are 
confined to the surface of a core. 

3. The effective surface concentration of adsorbed sub- 
units is increased if subunits are confined within less 
than a molecular diameter a by the core-subunit inter- 
action, so we take p^"^ = ^3(jv^'_„) exp (S'sd/fefi). This 
factor ensures that entropy losses due to subunit local- 
ization are not counted twice. 

4. Adsorption rate, fcad. The rate of subunit adsorption to 
core surfaces in our simulations can be calculated from 
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the Smoluchowski equation, with forces due to the ad- 
sorption potential explicitly included. Instead, we use 
the diffusion limited rate for the zero force case, but 
take fcad = 47r_Di?cut, where i?cut = 2.5a is the subunit- 
core center of mass distance at the point when the in- 
teraction force becomes nonzero. This choice was val- 
idated by comparing theory and simulation results for 
the time dependence of the number of adsorbed sub- 
units in the case of no assembly, i.e. Sb = 0. Since 
the excluded volume of adsorbed subunits is accounted 
for in Eq. IA2I we set the adsorption blocking factor 
4'n.m — 1 (see Eqs. \T2\ . This approach slightly over 
estimates the adsorption rate at moderate coverage, but 
yields the proper equilibrium surface density if there is 
no assembly. Because of the definition of the unit time 
to, the subunit diffusion constant is D = 1/48. 

We find that the agreement in most cases is surprisingly 
good, considering the degree of error in estimating Siot and /, 
and the approximation that there is only one reaction pathway. 

APPENDIX B 

In this appendix we summarize the application of the 
method of Scheutjens and Fleer [S] to polyelectrolytes. Our 
presentation closely follows that given in Bohmer et al. 16711 . 
except there is a different geometry, the Flory-Huggins pa- 
rameter for interactions between polymer segments is set to 
X = 0, and there are additional terms in the free energy equa- 
tions. 

We consider a lattice around an impenetrable sphere with 
radius Rc, contacting bulk solution in layer M + 1, with 
M large enough that the presence of the surface is negli- 
gible in layer z — M. All quantities are assumed uni- 
form within a layer z. Fixed numbers of polyelectrolytes and 
TEG molecules are grafted in layers Zp = Reap /I = 26 and 
^TEG = Rc/l + 1 ~ 17, respectively, with the lattice dimen- 
sion / — 0.35 nm. 

In the calculation, a segment of species i is subjected to a 
potential field Ui{z), which is given by (relative to the bulk 
solution) 



Ui{z) ~ u'{z) + ViC'^^z) 



(Bl) 



volume fraction for all species sums to unity in layer z. The 
second term on the right-hand side accounts for electrostatic 
interactions, with the valency Vi for a segment of species i, e 
the charge on an electron, and ^'(z) the electrostatic potential 
in layer z. The statistical weight of a segment for species i in 
layer z, relative to the bulk solution, is given by the segment 
weighting factor 



Gi{z) = e-xp{~Ui{z) / UbT) 



(B2) 



so that the volume fraction of a free monomer of type i is 
given by 0i(z) = (l3\Gi{z). 

The volume fractions for segments in chain molecules are 
calculated with a segment distribution function Gi{z,s\l), 
which gives the statistical weight of a chain conformation 
starting with a segment 1, located anywhere in the lattice, and 
ending in layer z after s — 1 steps. The segment distribution 
function is calculated from a recurrence relation: 

G,{z,l\l) = G,{z) 

G,{z, s\l) - G,(z)[A_iG,(z - 1, s - 1|1) 

+Ao(G,(z, s - 1|1) + AiG,(z + 1, s - 1|1)] (B3) 

with A_i, Ao, and Ai the fraction of contacts for a lattice site in 
layer z with respective layers z — l,z, and z + 1 (see Ref. 1 89]). 
The statistical weight Gi{z, s\r) of a chain conformation that 
starts at segment r and ends with segment s in layer z is cal- 
culated in the same way. The segment distribution functions 
for grafted molecules are modified to constrain segment 1 to 
begin in the grafted layer as described in Ref. i90l] . To ac- 
count for impenetrable surfaces, Ui{z) = oo for z < zjeg for 
all species, and Ui{z) = oo for z > Zp for polyelectrolyte and 
TEG molecules. 

The volume fraction of segment s in a chain of r segments 
is determined from the joint probability that a chain conforma- 
tion starts at segment 1 and ends at segment s in layer z, and 
a chain conformation starts at segment r and ends at segment 
s in layer z: 



s) = G,G,(z, s|l)G,(z, s|r)/G,(z) 



(B4) 



where the term u'{z) represents the hard-core interaction, 
which is the same for every species and ensures that the total 



where the probabilities are divided by Gi{z) to avoid double 
counting of segment s and Gi is a normalization constant. If 
molecules of species i are free to exchange with the bulk so- 
lution, Gi is determined by the fact that all segment weighting 
factors Gi = 1 in bulk so that 



G^ = <^\It, 



(B5) 
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If the total amount of molecules of species i is fixed (i.e. for a 
grafted brush) Ci is determined from 



a=0^/[nG,ir\l)] 



(B6) 



with the total amount of molecules of species i calculated 
from 

M 

0, (B7) 

and the chain weighting factor Gi{r\l) — 
^^j^ L(z)Gi(z, r|l), which gives the statistical weight 
for a chain of type i to be found anywhere in the lattice. The 
number of lattice sites L{z) depends on the layer z because 
of curvature and is calculated in Ref. ll89ll . 

The electrostatic potential. Electrostatics are accounted for 
with a multi-Stern-layer model in which all charges within a 
lattice layer are located on the plane at the center of the layer, 
with no charge outside of that plane. The plane charge density 
in layer z, a{z) is calculated from a sum over all species 



a{z) 



\ai(z)vie(j)i{z)/a^ 



(B8) 



where we set the cross-sectional area per lattice site — P 
with / the distance between layers. For weak acid groups the 
degree of dissociation is given by ll68l 



ai{z) 



(B9) 



\</'H30+(2)) +^»('/'H20(^)) 

with the dimensionless dissociation parameter Ki related to 
the dissociation constant by Ki ~ jci with the mo- 
larity of pure segments for species i. 

The electrostatic potential with respect to bulk solution at 
layer M is calculated through electroneutrality [67] and the 
potential in each layer is calculated from Gauss' law, for 
spherical layers it is 16811 



+ 1) 



1 



e(z + l)z(z + l/2) 



e(z)(z - l/2)z 

Z 

L{z')a{z') (BIO) 

z' = l 



with e{z) the volume fraction averaged dielectric permittivity 
for layer z. 

Free energies. The free energy is obtained from the canon- 
ical partition function for the lattice system 11661 h 
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m 



A = kuT 



^-ln(r,C,) 



}^ }^ L{z)c^,{z)u,{z) + JJel + i^diss (Bll) 



with the electrostatic energy for the case of a fixed surface 
charge given by 

M+l 

(B12) 



EL 



, M+l 

= - ^ L{z)a{z)^(z). 



2 = 1 



The final term in Eq. IB 11 1 was omitted in Ref. ll67ll . but is 
necessary to describe the free energy due to dissociation of 
acid groups 



M 



z—l a 



aa(2)log 



^H30+ 



(B13) 



aa{z)\ogaa{z) + (1 - aa{z))\og{l - aa{z)) 

where the sum extends over the weak acid species a. 
The Gibbs excess free energy A^x is obtained from 



A. 



A 



(B14) 



where the sum extends over the species i' that are in equilib- 
rium with the bulk solution. The Flory-Huggins equation for 
the chemical potential fii of a mixture of chain molecules is 
161 (withx = 0) 

For calculations with a fixed density of grafted polymers, the 
free energy includes additional terms, and the free energy /sf 
(per area at the interior surface of the capsid) is 



/sF({o-g})-^(Zp)Q.s 
fcfiT 



-4^+^L(Zg)[0-gl0gCrg 



+ (1 - Ug) log(l - ag) + agivg - l)] (B 1 6) 

where the sum is over the g £ {i} \ {i'} grafted species and 
Cg — Qgj {L{zg)rg) Is thc graftlug density (per lattice site) of 
species g. The first two terms in the sum account for fixed 
locations of grafted (or assembled in a capsid) segments ll92ll 
and the third adjusts the polymer reference state from the pure 
liquid to solvated polymer in infinite dilution (the last two 
terms of Eq. lBTSl with cj)^ = 0). With this choice of reference 
state, the excess free energy /ex as a function of the amount 
of adsorbed polyelectrolyte is 



/Ex(crp) = /sf(ctp) + A:BTCTplog(a^pp) 



(B17) 



where is the standard state volume and we have suppressed 
the dependence on ctjeg, which is constant for our calcula- 
tions. 
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We consider 9 molecular species, H2O, OH~, H3O , 
monovalent positive salt ions, monovalent negative salt ions, 
divalent positive salt ions, polyelectrolyte, neutral TEG, and 
acidic TEG end groups. Following Bohmer et al. 11671] the 



number of species is reduced to 8 by treating water as a 
weak electrolyte with valence —1 and dissociation constant 
pK^ = 14 - 21og(103iVA/^) ( Eq.lUfor 0w = 1, Avo- 
gadro's number, and P the volume of a lattice site). 
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